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The aim of this paper is the formulation of the finite element method in polar coordinates to solve transient heat 
conduction problems. It is hard to find in the literature a formulation of the finite element method (FEM) in polar 
or cylindrical coordinates for the solution of heat transfer problems. This document shows how to apply the most 
often used boundary conditions. The global equation system is solved by the Crank-Nicolson method. The pro- 
posed algorithm is verified in three numerical tests. In the first example, the obtained transient temperature dis- 
tribution is compared with the temperature obtained from the presented analytical solution. In the second numeri- 
cal example, the variable boundary condition is assumed. In the last numerical example the component with the 
shape different than cylindrical is used. All examples show that the introduction of the polar coordinate system 
gives better results than in the Cartesian coordinate system. The finite element method formulation in polar 
coordinates is valuable since it provides a higher accuracy of the calculations without compacting the mesh in 
cylindrical or similar to tubular components. The proposed method can be applied for circular elements such as 
boiler drums, outlet headers, flux tubes. This algorithm can be useful during the solution of inverse problems, 
which do not allow for high density grid. This method can calculate the temperature distribution in the bodies of 


different properties in the circumferential and the radial direction. The presented algorithm can be developed for 


other coordinate systems. The examples demonstrate a good accuracy and stability of the proposed method. 
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Introduction 


Solving of direct and inverse heat conduction prob- 
lems in cylindrical coordinate systems has received con- 
siderable interest, because of its extensive industrial ap- 
plicability in components, such as boiler drums, outlet 
headers, pipelines. A formulation of the control volume 
method in cylindrical coordinates to describe melting 
process is presented in [1]. Solution of the hyperbolic 
heat conduction problems in the cylindrical coordinate 
system by a hybrid Green’s function method can be 
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found in [2]. Formulation of the control volume method 
for the solution of an inverse heat conduction problem in 
the steam header cross section is shown in [3]. It is hard 
to find in the literature a formulation of the finite element 
method (FEM) in polar or cylindrical coordinates for the 
solution of heat transfer problems. Many finite element 
programs allow the geometry to be described in cylin- 
drical or spherical coordinate systems, but the elements 
in these codes are usually formulated in Cartesian coor- 
dinates. Exceptions are axisymmetric finite elements. 
Some studies have been done for development of solid 
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finite elements in spherical or cylindrical coordinates. 
Three-dimensional finite element analysis in cylindrical 
coordinates for nonlinear solid mechanics problems was 
presented in [4]. Full-vectorial FEM in a cylindrical 
coordinate system for loss analysis of photonic wire 
bends can be found in [5]. The aim of this paper is to: (a) 
describe the polar coordinate finite element formulation 
for transient heat conduction; (b) outline boundary 
condition formulations; and (c) show the effectiveness of 
the formulation for three selected problems. 


Formulation of the FEM 


Three-dimensional transient heat conduction analysis 
in many components of a power unit, such as: turbine hou- 
sings, pipelines, boiler drums, headers is often conducted 
in the chosen cylindrical cross-section. Two-dimensional 
temperature distribution can be described in the polar 
coordinate system as T(r, 0, ô. 

The transient heat conduction equation can be written 
in the following matrix form [6] 

op =-{g}" {G}+ay (1) 
where c, p are specific heat and density, {g} is the gra- 
dient column vector operator, {g}is the heat flux vector 


and gy is the heat generation rate per unit volume. The 
term fey {q} should be understood as V-{q}, where 


(V -) is the divergence operator. 
The gradient vector operator in the polar coordinate 
system (r, @) has the form 


A 
or 
— 2 
t=) a (2) 
r 00 
The heat flux vector q is defined by Fourier’s law 
(= -[Pligh? (3) 


where [D] is the conductivity matrix and in 
two-dimensional polar coordinates, it has the form 


e | (4) 


0 &(T) 
For isotropic material k,(T)=k,(T)=k(T). All ma- 


terial properties (c, p, k) are assumed as known functions 
of temperature. 

Transient heat conduction problems are initial-boun- 
dary problems for which one is required to assign appro- 
priate initial and boundary conditions. Initial conditions, 
also called Cauchy conditions, are temperature values of 
a body at its first moment tọ = 0 s. 


T(r,0.t)| o =i (r,8) (5) 


Three the most often used boundary condition of 1st, 


2nd and 3rd order can be assigned on the body’s edge 


T|, =7 (6) 
(OEDI M G 
(a t), =+(7-7,) 8 


å s 


Fig. 1 A diagram with different boundary conditions. 


Ta 


where 

{n} — unit outward normal vector to the boundary /; 

T, — temperature set on the body boundary 17, 

å, — heat flux set on the body boundary 7} , 

h — heat transfer coefficient set on the body boundary 
Th, 

Tm— temperature of a medium. 

Initial-boundary problem (1,5-8) was formulated for 
the whole region 2. Temperature can very in time and 
space and its distribution inside the finite element 2 ° is 
approximated by the following function 


T(r,0) = ÈT N, (é.7)={N}" {T°} (9) 


where n is the number of nodes in the element, T? — 
temperature in j-node and N(¢, 7) the shape function 
(interpolation function), {7°} is the element nodal 
temperature vector, {NV} is the shape function vector. For 
the four node quadrilateral element shape function vector 
is formulated as follow 
(l-—g)U—7)/4 
1+€)1-7)/4 
py =f C90- a 
(1+&)+7)/4 
(l—¢)U+)/4 
The time derivative of equation (9) has the form 


ThA. T(r,0) = LIN, (En)={n} {re} ay 


The temperature gradient vector VT inside the ele- 
ment 2° can be calculated from the following equation 

{g}7 =[B] {T°} (12) 

where [B] is the 2 x 4 matrix, which contains partial 


derivatives of the shape functions with respect to the 
global r and @ directions. 
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[B]={g}{y}" (13) 

Evaluating the [B] matrix from this equation can be 

problematic, because the shape function vector {N} is 

defined in the normalize local coordinate system (é, 77) 
presented in Fig. 2. 


2 
r 


Fig. 2 Four node quadrilateral finite element in global and 
normalize local coordinates. 


It is better to use the following equation 
[B]=[J]' [DER] (14) 
where the matrix [DER] contains derivatives of shape 
functions with respect to normalize local coordinates (€, 77) 


ON, ôN, ON, 

[DER] = a - (15) 
ôN, ôN, ON, 
an Ope F 


[J] is the transformation matrix between the global (r, 
0) and the normalize local (é, 17) coordinate system pre- 
sented in Fig. 2. 


a „20 
|3 06 

[J]= or a0 (16) 
ôn On 


where r can be calculated inside the element by mul- 
tiplying the shape function vector {N} and the vector of 
radial coordinates of element nodes 


r={Ny" (17) 


The partial derivatives of r, 0 , in respect to é 7 can 
be found by differentiation of coordinates expressed 
through shape functions and nodal coordinate values: 

or ON; ON; 


a0 
= ; = i @. 
ae ae ge Loe? an 
ar _ ON, 00_ AN, 
ee iy _ ig 
én a én i ôn a dn ' 


Using the principle of virtual temperatures the set of 
equations (1) and (5-8) can be reduced to the following 
set of equations [7] 
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[oe] e =f} fr) far] 


where [M* |= J pe{N}{N}'dédn (20) 
[xé] -J [8] [D] [4] gan (21) 

[xi] fay a (22) 

lh fate (23) 

(E= [anti (24) 

ji (= Jonna 05) 


If a convective boundary condition occurs at the edge 
1-4 (Fig.2), Eq. (22) has the form 


[Ki] [a-{ve=-p} 


2 2 
T || Or 00 
{N(é=-D} Hr dn 
on on 
er re are components of the transformation 
7 


on 
matrix presented in eq.(16) [12]. The matrix Ki | for 


(26) 


plane elements can be approximated as a diagonal matrix, 
with the diagonal terms defined as 


ary ( ao” 
h-IN(£=-1 pl gee ¢ 27 
Lee nyse) oe] 


ôn 
Equation (19) can be simplified as 
[m t irre 
A generalized Crank-Nicolson method, also known as 
0 method, can be applied to integrate numerically Eq.(28). 


Between the temperatures at the time t”"' and t” = nAt, n 
= 0, 1, ... the following dependence occurs 


{r° ye ={re\" + (1 -o){re\" +0{re | At (29) 


where 0 < 0 < 1. Eq.(29) is known as a generalized 
trapezoidal approximation [8]. Eq. (28) can be written for 


t” and t”™' 
[Ey ry Ey 
MIETET E o» 
The first set of equations should be multiplied on both 
sides by (1—0 ), while the second by@ . By adding the 


sides of equations (30) and (31) and after simple trans- 
formations one gets 


(28) 


(30) 
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(m J a[i = 
[m ]-(1-0)a [x ]]{ry + 
(1-0)Ar{ FY + oarre y” 
If 0 > 1/2, then the solution stability is ensured for the 
arbitrary time step At. However, time step At should be 
small due to the accuracy of temperature determination. 


In this paper, the unconditionally stable Crank-Nicolson 
method, was used assuming 0=1/2. 


(32) 


Analytical solution for one-dimensional transient 
temperature in a long pipe 


A long pipe with an inner radius r; and an outer radius 


ro has initial uniform temperature distribution of Ty. Next, 


the inner and outer surfaces of the tube are exposed to 
surrounding media at constant temperature T„ and T,, 
respectively. Heat transfer coefficients on the inner and 
outer surfaces are defined as h; and ho. 

The presented initial-boundary problem can be de- 
scribed by the equations (1,5,8). Assuming that q;=0, tem- 
perature independent material properties and 7=7(r, f), 
the set of equations (1,5,8) simplifies to 


oT 10 oT 
=k — 33 
Oo A m) on 
Tria (34) 
ot oa (T-7,) on =r; (35) 
or 
kanit) on r=r (36) 
or 


Applying the Laplace transform with respect to ¢ to the 
above equations, the solution can be written in the form [9] 


+ as 
hr; hor, 


ry ( o) f (tnr) 


"=n? +k u -G7 (he +k us) 


T,,h; -T,G,,h 


m i 


evn + (37) 


Li ACO 
2 i HP +k u -G3 (he +k) 
f’ -Tf (tns)ndn 
where u, are the positive roots of the equation 
[ro (4,75) + ke (aar) 
[Ao (Ant ) -k4p Yi (Ant )] 
[hoo (Hats) -ktn (tara) |: 
LAY (Han) +k MAY (tta; )]=0 


(38) 
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and f(z, r) and G, are given by 
S (tnr) =Y (Hats) + kttpYi (Hats) Jo ar) - 6) 
-| h;Jo (Hat) + kti (eats) Yo (Har) 
hYo ( k 
[a 0 (u, )+ HY, (44,7;) | (40) 


[haYo Hat) —KenYi (att) | 

Jo(ur) is the Bessel function of the first kind of the 
zero order. Yo(z) is the Bessel function of the second 
kind of the zero order, Jı(ur) is the Bessel function of the 


first kind of the first order. Yı(ur) is the Bessel function 
of the second kind of the first order. 


Numerical examples 


The outer and the inner radius of a long horizontal cy- 
lindrical header equals r= 177.8 [mm] and r= 127.8 
[mm] respectively. It is made of martensitic P91 steel. 
Thermal properties of the P91 steel are taken at the tem- 
perature of 100 [°C] and are: c= 486 [J/kgK], p= 7750 
[kg/m*], k=29 [W/mK]. An initial temperature of the 
header equals 7)=19.4 [°C]. Next, its inner and outer sur- 
faces are exposed to the medium 7,,=146.8 [°C] and 
T,=19.4 [°C] respectively. The heat transfer coefficients 
on the inner and outer surfaces are assumed as h; = 2000 
[W/mK] and h,=2 [W/mK] based on research done 
[10]. The transient temperature distribution is obtained 
from analytical solution presented in chapter 3. For the 
assumed data, the following six positive roots of Eq. (38) 
are obtained: 4 = 22.719, 4b = 77.174, 16 = 135.107, 44 = 
195.327, us = 256.626, L= 318.469. Temperature calcu- 


lated from Eq. (37) is presented in Fig. 4 and Table 1. Fig. 


4 shows temperature transients at the points P1 and P2. 
Location of these points is illustrated in Fig. 3. Exact 
values of temperature after 100 [s] of the heating opera- 
tion at selected five points through the wall thickness are 
presented in Table 1 as “T — exact”. 

Next, transient temperature distribution is calculated 
by FEM using the time step of A1 [s]. 

First, the header cross section is divided into four node 
quadrilateral finite elements formulated in Cartesian 
coordinates. The division into four elements across the 
wall thickness and four elements along the circumference 
is presented in Fig. 3. The similar division into six finite 
elements along the circumference and three across the 
wall thickness is often proposed during the solution of 
inverse heat conduction problems in cylindrical elements. 
Further increasing of the grid density may lead to 
instability of the inverse solution [11]. The results ob- 
tained are marked in Table 1 as “T - num. Cart.” 

Relative errors E between the exact and the calculated 
by FEM temperature values are calculated as 


pale 7 [Fer -Tren |100% 


ex 


(41) 
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Fig. 3 Finite element discretization in Cartesian and polar 
coordinates, location of the points P1, P2, P3, P4, PS. 
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Fig. 4 A comparison between the exact temperature transients 
at points P1, P5 and the temperature transients calcu- 
lated by FEM in Cartesian and polar coordinates. 


where Tx stands for the exact value and Trpm for the 
calculated temperature. The maximum error occurs on 
the outer surface at the node P5 and equals after 100 [s] 
about 11%. The relative error as the function of time on 
the inner and outer surfaces can be seen in Fig.5. 

A significant increase in accuracy can be achieved by 
using quadratic quadrilateral element with eight nodes 
without changing the coordinate system (T - num. Cart. 
nl.). However, if the number of linear equations is to be 
maintained constant, the number of elements should be 
smaller. This results in the increase of the error close to 
the inner surface at the node P1. 

The formulation of FEM in the polar coordinate sys- 
tem allows for the use of a linear four node quadrilateral 
finite element and allow to achieve better accuracy at all 
points (T - num. pol.). 
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Fig. 5 The relative error E between the exact temperature 
transients at points P1, P5 and the temperature tran- 
sients calculated by FEM in Cartesian and polar coor- 
dinates. 


To receive better accuracy in Cartesian coordinates 
than in polar coordinates the mesh should be denser. The 
division into 16 elements across the wall thickness and 
16 elements along the circumference is marked in Table 1 
as “T - num. Cart. fine mesh ”. This discretization will be 
used in the next numerical example as an exact solution. 

In the second numerical example the variable boun- 
dary condition is assumed to verify the proposed boun- 
dary condition formulation in the polar coordinate system. 
The 3rd kind boundary condition, presented in Eq.(8), is 
the most difficult to apply because it requires the defini- 
tion of the matrix in Eq.(22) and the vector in Eq. (25). 

The medium temperature inside the header is assumed 
to be constant Tm =146.8 [°C[. The heat transfer coeffi- 
cient on the inner surface is also constant in time but va- 
ries depending on the angle with the following function 

h(0)=1650-1350-cos(8) (42) 

It means that A at the highest point of the header 
equals h(@=180°) = 3000 W/m?°K and at the lowest point 
h( @=0°) = 300W/m’K. 

The temperature obtained after 100 [s] of the heating 
operation at the same five points as in the first example is 
presented in Table 2. The results achieved by four node 
quadrilateral finite element in the polar coordinate sys- 
tem (T - num. pol.) are the most similar to the exact solu- 
tion. 

The last numerical example is intended to show that 
the proposed formulation of FEM in the polar coordinate 
system can also be applied to elements of shapes other 
than cylindrical. The geometry of the analyzed compo- 
nent in Cartesian and polar coordinates is presented in 
Fig. 6. The proposed shape is similar to the geometry of 
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the flux tube used during the identification of local heat 
flux to membrane water-walls in steam boilers [13]. The 
component outer surface is the same as in the previous 
examples, but the inner surface was moved 10 mm up. 
The assumed boundary conditions are the same as in the 
second example. 

The results obtained in the polar coordinate system (T 
- num. pol.) are the most similar to the exact solution. 
Only on the outer surface at the point P5 the solution 
obtained by quadratic quadrilateral element is better. 


Conclusions 


The finite element formulation in the polar coordinate 
system for the solution of transient heat conduction 
problems has been presented. The proposed algorithm 
has been verified in three numerical tests. In the first 
example, the obtained transient temperature distribution 
is compared with the temperature obtained from the pre- 
sented analytical solution. The solution improvement 
caused by the introduction of the polar coordinate system 


Table 1 The exact and the calculated temperature after 100 [s] of the heating operation at five points through the wall 


thickness of cylindrical component 


Temperature location Pl P2 P3 P4 P5 
T - exact [°C] 110.32 82.27 61.81 49.58 45.56 
T - num. Cart. [°C] 110.46 84.66 65.73 54.35 50.61 
T - num. Cart. nl. [°C] 107.87 80.64 60.97 49.06 45.33 
T - num. Cart. fine mesh [°C] 110.27 82.29 61.89 49.69 45.69 
T - num. pol. [°C] 110.21 82.08 61.50 49.16 45.13 
Error E num. Cart. [%] 0.13 2.90 6.35 9.62 11.09 
Error E num. Cart. nl. [%] 2.22 1.98 1.37 1.05 0.51 
Error E num. Cart. Fine mesh [%] 0.05 0.03 0.13 0.23 0.29 
Error E num. pol. [%] 0.10 0.23 0.51 0.84 0.94 


Table 2 The exact and the calculated temperature after 100 [s] of the heating operation at five points through the wall 
thickness of cylindrical component with variable boundary condition. 


Temperature location Pl P2 P3 P4 P5 
T - exact [°C] 104.23 77.51 58.27 46.87 43.15 
T - num. Cart. [°C] 105.11 80.30 62.31 51.56 48.07 
T - num. Cart. nl. [°C] 100.70 75.07 57.32 46.11 42.88 
T - num. pol. [°C] 105.01 77.98 58.42 46.78 42.99 
Error E num. Cart. [%] 0.84 3.60 6.92 10.00 11.41 
Error E num. Cart. nl. [%] 3.39 3.14 1.64 1.62 0.62 
Error E num. pol. [%] 0.75 0.61 0.26 0.19 0.36 


Table 3 The exact and the calculated temperature after 100 [s] of the heating operation at five points through the wall 
thickness of non-cylindrical component with variable boundary condition. 


Temperature location P1 P2 P3 P4 P5 

T - exact [°C] 104.28 77.12 57.67 46.14 42.20 
T - num. Cart. [°C] 105.25 79.71 61.18 50.00 46.05 
T - num. Cart. nl. [°C] 100.71 74.84 56.96 45.51 42.02 
T - num. pol. [°C] 105.12 77.44 57.46 45.51 41.44 
Error E num. Cart. [%] 0.93 3.36 6.08 8.35 9.13 
Error E num. Cart. nl. [%] 3.42 2.96 1.23 1.37 0.42 
Error E num. pol. [%] 0.80 0.42 0.38 1.37 1.80 


Fig. 6 Finite element discretization of the component in Car- 
tesian and polar coordinates. 


has been shown in time and the time of 100 [s] has been 
chosen for the comparison of temperature values ob- 
tained from different methods. 

The maximum relative error of the temperature distri- 
bution calculated by FEM in Cartesian coordinates oc- 
curs on the outer surface and equals about 11%. The pro- 
posed algorithm of FEM in the polar coordinate system 
allows to reduce the maximum relative error to about 1%. 
It has been presented that the introduction of the polar 
coordinate system gives better results than the use of 
higher-order elements, when the tasks with a similar 
number of equations are compared. 

In the second example, the variable boundary condi- 
tion has been assumed and the reduction of errors is even 
greater, because of the higher circumferential heat flow 
than in the previous example. Polar coordinates allow to 
model this phenomenon more correctly, than Cartesian 
coordinates. The last numerical example is the proof that 
the proposed algorithm can be applied not only to the 
cylindrical components. 

FEM formulation in polar coordinates is valuable 
since it gives a higher accuracy of the calculations with- 
out compacting the mesh in solids of circular shapes. 
This method can be useful during the solution of inverse 
problems, which do not allow for high density grid. This 
method can calculate the temperature distribution in the 
bodies of different properties in the circumferential and 
the radial directions. The presented algorithm can be 
developed for other coordinate systems. The applicability 
of the proposed method has been demonstrated by a way 
of comparison with the results obtained from the analyti- 
cal and numerical solution. The examples showed de- 
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monstrate a good accuracy and stability of the proposed 
method. 
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